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Abstract 

We obtain precise values for the fugacities of vortices in the 2-d 
planar rotor model from Monte Carlo simulations in the sector with no 
vortices. The bare spinwave stiffness is also calculated and shown to 
have significant anharmonicity. Using these as inputs in the KT recur- 
sion relations, we predict the temperature Tc = 0.925, using linearised 
equations, and Tc = 0.899 it .005 using next higher order corrections, 
at which vortex unbinding commences in the unconstrained system. 
The latter value, being in excellent agreement with all recent determi- 
nations of Tc, demonstrates that our method 1) constitutes a stringent 
measure of the relevance of higher order terms in KT theory and 2) 
can be used to obtain transition temperatures in similar systems with 
modest computational effort. 



* on leave from: Material Science Division, Indira Gandhi Centre for Atomic Research, 
Kalpakkam 603102, India 



1 



Introduction: The phase behaviour of isotropic magnets and related sys- 
tems in two dimensions is a particular challenge since famous theorems ex- 
clude long-range order but nevertheless phase transitions occur in models 
such as the two- dimensional XY ferro magnet [0 or the planar rotor model 
(S'f = cos0j, = sin0j, where Si is a two component spin at the site i with 
unit magnitude and orientation < 0j < 27r ), 

m = -^i: cos(0,-0,), (1) 

<ij> 

the sum extends once over all nearest neighbor pairs of the (square) lattice, 
and T is the reduced temperature (the Boltzmann constant ks = /?~^T 
is taken to be 1 throughout). Originally |0| it was proposed that a critical 
temperature Tc occurs where the correlation length ^ describing the decay 
of the correlation function g{r) =< S(0)S(r) > with distance r, and the 
susceptiblity x = J2r < S(0) ■ S(r) > /T diverge according to power laws, 

^^t-',xr-t-\t = T/T,-l (2) 

z/, 7 being the usual critical exponents. However, Kosterlitz and Thouless 
(KT) 0, H developed a completely different scenario, based on the unbinding 
of vortex- antivortex pairs, yielding an essential singularity, 

Ine = In^o + ^r^ (3) 

where b are nonuniversal constants, while an approximate renormalization 
group treatment]^, ^ predicted that the exponents u,ri take the universal 
values, 

u = 1/2, (4) 
r] = 1/4. 

For T < Tc, spin- wave theory remains essentially valid and g{r) ~ r^^^"^) 
where ti{T) increases smoothly with increasing temperature from ?7(T = 0) = 
upto ?7(T'c) = ?7 = 1/4. The spinwave stiffness K{T) (in which we have ab- 
sorbed a factor of 1/T as in Eq.(l)) smoothly decreases and also involves 
a universal ratio at Tc,K{Tc) = 2/7r[^, ^. A related critical behaviour 
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is predicted for the superfluid- normal fluid transition of helium in two- 
dimensions]^], for the roughening transition of interfaces 0, transitions of 
adsorbed layers on surfaces to modulated structures incommensurate with 
the substrate periodicity, etc. Thus, this problem has found widespread 
interest [0]. A particularly interesting — - but also still controversial — - ex- 
tension deals with two- dimensional melting p!0[]. 

However, both the physical mechanism for the vortex- antivortex pair 
dissociation at Tc and the resulting predictions have been questioned many 
times (eg.[0, H). The theory|, |, |, |15|, |16], 0, H involves problem- 

atic assumptions such as the decoupling of vortex and spinwave excitations; 
and numerical analyses |19|, |2^, 



|, 0, ^ often are not fully 
convincing although usually the KT theory is favoured. Monte Carlo studies 
are difficult since ^ increases so strongly as t gets small (eg. ^ > 40 lat- 
tice spacings for t < .1), and so it is questionable whether the asymptotic 
critical region is reached. Even studies for very large lattices (1200 x 1200) 
still reveal problems with Eq.(4)^3|], and the simulation data can well be 
fitted to Eq.(2) if (albeit rather large) corrections to scaling are taken into 
account [^. The most recent analyses, in fact, point towards the need of 
considering logarithmic corrections p6| . 



In the present paper we hence follow a different strategy for a Monte Carlo 
test of the KT renormalization approach, avoiding the brute force methods 

j, n, n, n, g 



of Refs.O, 19, 20, 21 



Namely we test the KT scenario 
by obtaining the proper input parameters for the renormalization group flow 
equations, which then are solved numerically. In this way a stringent consis- 
tency test is possible which is different from all previous approaches to the 
problem. 



Monte Carlo estimation of input parameters to the KT theory The 

KT theory p, |[ |^, |^, |^ can be cast in the framework of a two parameter renor- 
malization flow for the spinwave stiffness K{1) and the fugacity of vortices 
y{l), where / is related to the considered length scale as / = ln(r/a), where a 
is the lattice spacing. These flow equations in terms of the scaled variables 
X = {2 — 7iK) and y' = Airy and upto next to leading order|l7] are, 

-ji = y -y-^ (5) 
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dy' , 5 ,o 

Using a linearised version of these equations (i.e. keeping only the first terms 
on the right hand side of these equations) and using the approximate initial 
conditions y{l = 0) ~ exp(— 10.2/2T) and K{1 = 0) = 1/T - a result from 
harmonic spin wave theory strictly valid at T — > 0, Kosterlitz[Q found that 
a non- trivial fixed point K{1 = oo) = 2/11, y{l = 00) = 0) exists (cf. Eq. (3) 
above) but the resulting estimate for Tc ~ 1.35 is rather different from the 
current best estimates Tc = 0.895 ± .005 |2^, |2^, Does this discrepancy 
mean that the KT scenario does not work? 

Such a conclusion would clearly be premature, however, because the above 
assumption implies that even at T = Tc one can still take the unrenormalized 
zero temperature value of the spin wave stiffness as a starting value for the 
recursion, Eq.(5). To test this assumption we have obtained K (and y) 
from Monte Carlo simulations. Two sets of simulations are carried out. The 
first set uses the full Hamiltonian, Eq.(l), while the second set uses the 



constraint that neither vortices nor antivortices can form 112 7||. Note that an 



elementary plaquette of the square lattice contains a vortex or an antivortex, 
if the angles (pi of the spins 1, 2, 3, 4 at the corners of the plaquette (labelled 
anticlockwise) satisfy the condition J2i=i ^4>i = ±27r, where we have defined 
A0i = — 0j,05 = 01. If there are only spin wave excitations in the 
system, X]f=i A0j = for all plaquettes. Hence the no- vorticity constraint in 
the Monte Carlo elementary step (which involves an attempt to replace (pi 
by a randomly chosen (p[, with < 0^ < 2tt) considers whether J2i=i ^(pi = 
is still true for this trial configuration for all the four adjoining plaquettes 
to which the site i belongs. If the constraint is not true, the trial move 
is automatically rejected. Note that we always start the simulation from a 
vortex free initial fully aligned state {cos (pi = 1 for all i). 

Fig.l gives a plot of the inverse stiffness constant = 47rln < > 
/InN, where < >= Txn- Note that while diverges in the uncon- 
strained system at T^, it stays finite in the constrained system and finite size 
effects are negligible even at Tc since the constrained system is not at a criti- 
cal point there. Therefore K can be obtained very precisely - the constrained 
system was equilibrated using 2 x 10^ Monte Carlo Steps (MCS) per site and 
a further averaging over 3 x 10^ MCS was sufficient to obtain high quality 
data. Fig. 1 shows that indeed the harmonic theory result for K^^{= T) is 
poor near Tc. 
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Next we wish to estimate y{l = 0) from simulations as accurately as 
possible. This was done in two ways: (i) the concentration of vortex pair 
excitations was measured in the unconstrained simulations of Eq.(l), (ii) the 
rejection rate p of Monte Carlo moves that were rejected was measured in the 
constrained simulations (Fig. 2.). The chemical potential n of the vortices 
could be obtained from the slope of — ln?7,„ or —Inp as a function of T~^. 
We see that again our data for p in the constrained simulations were of much 
superior quality because of the absence of a phase transition. Our estimate 
for y{l = 0) = exp(— yu/T) can now be used together with our value for 
K~^{1 = 0){— T + from numerical fits to the data) in the recursion 
relations Eq.(5) to obtain the renormalized rigidity modulus Kr and hence 
T 

The recursion relations Eq.(5) are solved numerically to obtain the renor- 
malized modulus and fugacity. Using only the linearized equations (which 
can be solved analytically) we obtain = 0.925 which is considerably closer 
to the experimental value than the KT estimate of 1.35 but there is still 
a significant discrepancy. However, this discrepancy vanishes when the full 
equations incorporating leading order correction terms (which do not affect 
universal behaviour and hence usually omitted) are used. Taking leading or- 
der correction terms into account, we obtain Tc = 0.899 ± 0.005 in excellent 
agreement with the brute force simulations! Fig. 3 presents the resulting flow 
diagram, which displays the importance of these correction terms directly. 

In addition, our results offer a simple way of calculating the nonuniversal 
critical amplitude b (note that Fig. 3 implies that Eq.(4) holds, of course). 
Critical amplitudes are usually notoriously difficult quantities to estimate 
directly from simulations. Following Ref.ljl^ we define yo as the intercept 
of the flows for T > Tc with the y'— axis (see Fig. 3). The leading order 
behaviour of the correlation length ^ ~ i^/yo- Using the fact that y^ has 
an expansion ?/q = J^iO'it^ with t = (T — Tc)/Tc (see Fig. 4), we get b = 
Ti I \fa^ = 1.534 ± .002 which is in excellent agreement with the estimate 
b = 1.585(9) obtained by 01ssonp4| by directly fitting the behaviour of the 
dielectric function e(t) = K/Kr. 



Discussion and conclusions Our analysis demonstrates that the planar 
rotor model is fully consistent with the KT theory, but for a quantitatively 
accurate description of the transition (especially for nonuniversal quantities) 
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it is indispensible that the higher order terms in the recursion relations Eq. (5) 
are taken into account. This finding imphes that far away from Tc significant 
corrections to Eq.(3) are expected — this offers an explanation why the direct 
simulations [jl2|, |21], ^ ^ have difficulties in extracting the correct 

critical behaviour unambiguously. In contrast, our method yields critical 
properties with comparatively modest computational effort. We expect that 
analogous methods can be applied to other models that are expected to 
show KT transitions; in fact, we are currently undertaking an extension of 
our approach to the controversial issue of two- dimensional melting. 
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FIGURE CAPTIONS 



Fig 1. The inverse spinwave stiffness K^^{T) plotted vs. temperature, for 
two types of simulations: (i) The unconstrained system of L x L lattices with 
L — 100 (O); (ii) The constrained (vortex free) system for L = 30(x),60(+) 
and lOO(n), respectively. The dashed straight line represents the result of 
harmonic theory {K"^ — T) while the full curve is a fit to the form T + aT^ 
with a = .50± .01. 

Fig. 2. Plot of — In(n^) and — ln(p) vs. the inverse temperature. The vortex 
concentration n„ (O) is calculated in the unconstrained simulation of Eq.(l) 
for L X L lattices (L = 100). The rejection ratio p in the constrained simula- 
tion was calculated for L = 60(+) and L = lOO(n) respectively. The dotted 
line is a fit to the latter data (near the transition temperature) yielding 
2/i = c = 6.55 ± .03. 

Fig. 3. Flows of a; = (2 — tiK) and y = exp{—fi/T) under the action of the 
renormalisation group starting from a set of initial conditions (o) obtained 
from our simulations of the XY model. The dotted lines {y — ±x) are 
the separatrix for the linearized flow equations valid for flows near the flxed 
point X = 0,y = 0, the thick lines are the actual separatrix for the nonlinear 
equations Eq. (5). Note that these curves separate flows that terminate on 
the critical line x < 0,y = (ordered phase) from flows towards y — > oo 
(disordered phase). Arrows give the direction of the flow. 

Fig. 4. Square of the intercept yo, of the flows (Fig. 3.) with the y'— axis as 
a function of t = (T — T^j/Tc for T > T^.. The data points (O) are fltted to 
a straight hue y' = ait. The critical amplitude b for the correlation length 
^ ~ exp(6rV2), is given by 6 = = 1.534 ± .002. 
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Figure 1: 

(Sengupta, Nielaba , Binder; Euro. Pliys. Lett.) 
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Figure 2: 

(Sengupta, Nielaba , Binder; Euro. Phys. Lett.) 
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Figure 3: 

(Sengupta, Nielaba , Binder; Euro. Phys. Lett.) 



11 



0.01 




12 



